Introduction to geocomputation with R


shorturl.at/jCX09

Duc-Quang Nguyen

2022/05/19

Simple features: sf

Supports common vector geometry types

Simple feature types fully supported by sf.

Source: Geocomputation with R by Robin Lovelace

sf can represent all common vector geometry types (raster data classes are not supported by sf): points, lines, polygons and their respective ‘multi’ versions (which group together features of the same type into a single feature). sf also supports geometry collections, which can contain multiple geometry types in a single object.

Source: Geocomputation with R by Robin Lovelace

Coordinate reference system (CRS)

How the spatial elements of the data relate to the surface of the Earth. That coordinate systems are a key component of geographic objects

Knowing which CRS your data is in, and whether it is in geographic (lon/lat) or projected (typically meters), is important and has consequences for how R handles spatial and geometry operations

Source: Geocomputation with R by Robin Lovelace

Carte choroplèthe par municipalité

library(tidyverse)
library(sf)
library(swissdd) # Données historiques et en temps réel des votations fédérales 
library(ggiraph) #ggiraph_0.8.2

#require("revealjs")

# install.packages(c("tidyverse", "sf", "swissdd", "reveal.js", "ggiraph"))

# devtools::install_github('davidgohel/ggiraph')

Load Swiss geo data

Data from Office fédéral de la statistique - Limites communales généralisées: géodonnées

muni <- st_read("shp/g2g22.shp")
## Reading layer `g2g22' from data source 
##   `/Users/duc-q.nguyen/Documents/ddj/HEIGVD/R/shp/g2g22.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2148 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2485424 ymin: 1075268 xmax: 2833837 ymax: 1295937
## Projected CRS: CH1903+ / LV95

Anatomy of a sf object

str(muni)
## Classes 'sf' and 'data.frame':   2148 obs. of  20 variables:
##  $ GMDHISTID: num  13256 11742 11801 11992 12249 ...
##  $ GMDNR    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ GMDNAME  : chr  "Aeugst am Albis" "Affoltern am Albis" "Bonstetten" "Hausen am Albis" ...
##  $ BZHISTID : num  10053 10053 10053 10053 10053 ...
##  $ BZNR     : num  101 101 101 101 101 101 101 101 101 101 ...
##  $ KTNR     : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ GRNR     : num  4 4 4 4 4 4 4 4 4 4 ...
##  $ AREA_HA  : num  791 1059 743 1360 653 ...
##  $ E_MIN    : num  2678111 2673826 2675745 2680394 2674866 ...
##  $ E_MAX    : num  2681155 2678584 2679642 2686463 2678972 ...
##  $ N_MIN    : num  1234561 1235208 1238997 1230185 1238053 ...
##  $ N_MAX    : num  1238544 1239338 1243159 1236411 1241034 ...
##  $ E_CNTR   : num  2679300 2676800 2677800 2682900 2676400 ...
##  $ N_CNTR   : num  1235700 1236800 1241000 1233100 1239000 ...
##  $ Z_MIN    : num  533 440 502 524 475 463 412 386 432 383 ...
##  $ Z_MAX    : num  887 749 715 915 745 621 504 476 658 520 ...
##  $ Z_AVG    : num  685 529 581 673 560 557 454 414 506 440 ...
##  $ Z_MED    : num  673 502 583 653 543 559 454 410 493 439 ...
##  $ Z_CNTR   : num  700 490 544 610 503 573 431 400 471 453 ...
##  $ geometry :sfc_MULTIPOLYGON of length 2148; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:24, 1:2] 2678219 2678560 2678206 2678196 2678259 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:19] "GMDHISTID" "GMDNR" "GMDNAME" "BZHISTID" ...

Quick sf plots

# Quickly plot the geodata by selecting only the geometry
muni %>% select() %>% plot
# This will color the map by the features
# plot(muni)

# ggpplot2 way to visualize geo data
ggplot(muni) +
  geom_sf() 

Basic choropleth of Swiss municipalities by their canton

ggplot(muni) +
  geom_sf(aes(fill = KTNR))

Improved choropleth Swiss of municipalities by their canton

ggplot(muni) +
  # Hide the boundaries' line, discretize cantons
  geom_sf(aes(fill = as.factor(KTNR)), colour = NA) +
  # Hide colour scale
  scale_fill_viridis_d(guide = 'none') +
  # Remove graticule
  theme_void()

Choropleth map by municipality of the federal ballot «Lex Netflix» 15.05.2022

Données historiques et en temps réel des votations fédérales & cantonales suisses:
swissdd

id vote «Lex Netflix»: 6550

Get the votation results by municipality

muni_2022_05_15_results <- swissdd::get_nationalvotes(
  geolevel = "municipality", from_date = "2022-05-15", 
  to_date =  "2022-05-15", language = "FR")

# str(muni_2022_05_15_results)
#  muni_2022_05_15_results$name %>% unique()

ln_df <- muni_2022_05_15_results %>% 
  filter(id == 6550) %>% 
  select(canton_name:mun_name, jaStimmenInProzent, gueltigeStimmen)

Bind to geodata and plot

ln_df <- left_join(
  ln_df %>% mutate(mun_id = as.numeric(mun_id)), 
  muni %>% select(GMDNR:GRNR), 
  by = c("mun_id" = "GMDNR")
  # To ensure the join returns a sf object
) %>% st_as_sf()

ln_map <- ggplot(ln_df) +
  geom_sf(aes(fill = jaStimmenInProzent), colour = NA) +
  scale_fill_viridis_c(option = "E", limits = c(0, 100)) +
  theme_void() +
  theme(legend.position = "top") +
  labs(subtitle = "% de oui à la votation fédérale «Lex Netflix»")
ln_map 

Polish the map with cantons

# Recreate the cantonal boundaries with a spatial join

cant <- muni %>% 
  group_by(KTNR) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  ungroup() 

plot(cant)

ln_map +
  geom_sf(data = cant, size = 0.2, colour = "white", fill = NA)

Interactive R charts with HTMLwidgets

Bring the best of JavaScript data visualization to R

  • With a line or two of R code, use javaScript visualization libraries at the R console, just like plots
  • HTML widgets can be used at the R console as well as embedded in R Markdown reports and Shiny web applications

htmlwidgets gallery

Some htmlwidgets:

Interactive map with ggiraph

ggiraph makes ‘ggplot’ graphics interactive

Additional aesthetics:

  • tooltip: tooltips to be displayed when mouse is over elements.
  • data_id: id to be associated with elements (used for hover and click actions)

Call function girafe with the ggplot object so that the graphic is translated as a web interactive graphics. ggiraph

# simplify geodata to make it a bit quicker
ln_df <- ln_df %>% 
   st_simplify(preserveTopology = TRUE, dTolerance = 100)

ln_mapi <- ggplot(ln_df) +
  geom_sf_interactive(aes(
    fill = jaStimmenInProzent, 
    tooltip = jaStimmenInProzent,
    data_id = mun_id), 
    color = NA) +
  scale_fill_viridis_c(option = "E", limits = c(0, 100)) +
  theme_void() +
  theme(legend.position = "top")
girafe(ggobj = ln_mapi)

Better tooltip

# create more informative HTML tooltip
ln_df <- ln_df %>% 
  mutate(
    tp = str_c("<h3>", mun_name, "<h3>",
               "<i>", canton_name, "</i><br>",
               "<b>",
               prettyNum(jaStimmenInProzent, 
                         decimal.mark = ",", digits = 3),
               "% oui</b>"
               ) %>% 
      # ggiraph hack to ensure there's no single quote
      str_replace_all("'", "&#39;")
  )

ln_mapi <- ggplot(ln_df) +
  geom_sf_interactive(
    aes(fill = jaStimmenInProzent, tooltip = tp, data_id = mun_id), 
    color = NA) +
  scale_fill_viridis_c(option = "E", limits = c(0, 100)) +
  theme_void() +
  theme(legend.position = "top") +
  geom_sf(data = cant, size = 0.2, colour = "white", fill = NA)
girafe(ggobj = ln_mapi)

sf cheatsheet

Exercise 1

Make a choropleth map of the cantonal results of last Sunday federal ballot «Lex Netflix» (2022-05-15) with ggplot2

Make it static first and interactive with ggiraph

Help: * The “shp” folder contains the swiss cantonal boundaries “shp/g2k22.shp” * You should check the doc of: swissdd::get_nationalvotes